RNAseq Differential Isoform Report

Differential expression for Cut vs Uncut

This is the differential expression analysis report for the comparison of gene expression in Cut versus the gene expression in Uncut.

Data is processed from transcript quantifications derived from Salmon.

Files are normalized and compared between conditions using DESeq2. For more information on the analysis you can review the code in this document and visit our training material at https://rockefelleruniversity.github.io

outdir <- file.path(basdir,
                    paste0(conditionA,
                                "_vs_",
                                conditionB,
                           "_Reports")
                    )
dir.create(outdir,showWarnings = FALSE,recursive = TRUE)


files <- dir(countDir,
             recursive = TRUE,
             pattern="quant.sf",
             full.names = TRUE)
tx2gene <- read.delim(tx2geneFile,sep=",")
tximported <- tximport(files,type = "salmon",tx2gene = tx2gene,txOut = FALSE)
colData <- data.frame(Sample_Unique=gsub("_salmon","",basename(dirname(files))),
                      Group=gsub("_salmon","",basename(dirname(files))) %>% gsub("_.*","",.),
                      salmonCounts=files) %>%
  column_to_rownames("Sample_Unique")


colnames(tximported$counts) <- gsub("_salmon","",basename(dirname(files)))
colnames(tximported$abundance) <- gsub("_salmon","",basename(dirname(files)))

The Input files

The input files from Salmon are shown in Table 1 alongside the sample names and any group information.

inputTable <- data.frame(files=file.path(basename(dirname(files)),
                         basename(files)),
                         colData) %>% set_rownames(gsub("_salmon","",basename(dirname(files))))
inputTableforDE <- data.frame(files=files,
                         colData) %>% set_rownames(gsub("_salmon","",basename(dirname(files))))
DT::datatable(inputTable)

Additional files required for the analysis of differential isoform abundance are shown in Table 2.

DT::datatable(t(data.frame(TranscriptToGene=tx2geneFile,
GTF=pathToGTF,
GO=GOsets)) %>% as.data.frame %>% rownames_to_column %>% set_colnames(c("FileType","FileName")) %>% mutate(FileName=basename(FileName)))

Comparison of Gene expression

Data is subset to the conditions of interest - Cut and Uncut and imported into DESeq2.

Contents of differential expression report directory

The output files

The report directory contains output from a differential analysis including -

• The Reports - PDF, HTML and interactive HTML reports.

• Differential Expression tables - Tables of significance and extent of change for Cut minus Uncut.

  • Excel - Cut_vs_Uncut.xlsx
  • CSV - Cut_vs_Uncut.csv
  • In HTML pages (selfcontained) - Cut_vs_Uncut.html

• Figures - The Figures directory contains all images seen in reports.

The reports

Code for analysis

All code used in this report can be found by clicking code buttons on side of reports.

  dds <- DESeqDataSetFromTximport(tximported,colData = colData ,design = ~Group)
  dds <- dds[,colData(dds)$Group %in% c(conditionA,conditionB)]
  colData(dds)$Group <- factor(as.character(colData(dds)$Group))
  dds <- DESeq(dds)
  normExprs <- rlog(dds)
DEtable <- results(dds,contrast=c("Group",conditionA,conditionB)) %>% 
                as.data.frame %>% 
                rownames_to_column %>% 
                filter(!is.na(padj)) %>% 
                arrange(pvalue) %>% 
                arrange(desc(abs(stat)))

# save(dds,file="~/Downloads/Heintz/dds.RData")

normCountToExport <- sapply(levels(colData(dds)$Group),function(x)rowMeans(counts(dds,normalized=TRUE)[,colData(dds)$Group %in% x,drop=FALSE]))
colnames(normCountToExport) <- paste0(colnames(normCountToExport),"_NormCounts")
normCountToExport <- data.frame(IDs=rownames(dds),
                                normCountToExport)

normCountToExport2 <- lapply(levels(colData(dds)$Group),function(x)counts(dds,normalized=TRUE)[,colData(dds)$Group %in% x,drop=FALSE]) %>% do.call(cbind,.)
colnames(normCountToExport2) <- paste0(colnames(normCountToExport2),"_NormCounts")
normCountToExport2 <- data.frame(IDs=rownames(dds),
                                normCountToExport2)

# normCountToExport <- merge(annoS,normCountToExport,by=1,all=TRUE)
# write.table(normCountToExport,file=file.path(outdir,paste0(ComparisonName,"_CountsPerGroup.csv")),sep=",",row.names = FALSE)
DEtable <- merge(DEtable,normCountToExport,by=1,all.x=TRUE,all.y=FALSE) %>% 
           merge(normCountToExport2,by=1,all.x=TRUE,all.y=FALSE) %>% 
            arrange(desc(abs(stat)))
require(org.Dr.eg.db)
DEtableanno <- AnnotationDbi::select(x = org.Dr.eg.db,keys = as.character(DEtable[,1]),keytype = "ENSEMBL",columns = c("ENSEMBL","SYMBOL"))
DEtable <- merge(DEtableanno,DEtable,by=1,all.x=TRUE,all.y=FALSE)  %>% 
            arrange(desc(abs(stat)))
colnames(DEtable)[1] <- "rowname"

Sample heatmap

The Heatmap in Figure 1 shows an overview of the correlation between gene expression profiles for the samples under investigation. Samples with greater similarity in gene expression profiles will be clustered closer together in the heatmap and have smaller values for between samples distances (indicated by darker blue).

require(pheatmap)
sampleDists <- dist(t(assay(normExprs)))

sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- colData(normExprs)$Group
colnames(sampleDistMatrix) <- colnames(normExprs)
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
# png(paste0(ComparisonName,"_sampleDist.png"))
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         annotation_col=as.data.frame(colData(normExprs)) %>% 
           dplyr::select(!contains("sizeFactor")) %>% dplyr::select(Group),
         col=colors) 

# dev.off()
graphics.off()

Principal Component Analysis - PCA

The Principal Component Analysis in Figure 2 displays the separation of samples along the major sources of variation in the data.

p <- plotPCA(normExprs,ntop=500,intgroup="Group")+theme_bw()+ggtitle(paste0("PCA for top500 most variable genes in ",conditionA," and ",conditionB," samples"))

Differential expression analysis

In this section we compare the log2 expression of genes in Cut minus log2 expression in Uncut using DESeq2.

Table 2 shows the results of this comparison with columns

• IDs - Entrez IDs for genes
• Symbols - Symbols for genes
• baseMean - Average normalized counts for gene across all samples.
• log2FoldChange - Log2 mean difference in expression for Cut minus Uncut.
• lfcSE - Standard error of log2FoldChange estimate
• stat - Wald statistic for log2FoldChange and lfcSE.
• pvalue - Significance of log2FoldChange for Cut minus Uncut.
• padj - Adjusted pvalue for log2FoldChange for Cut minus Uncut to account for the multiple comparisons performed in testing differential expression across all genes).
• Cut_NormCounts - Normalised counts in Cut.
• Uncut_NormCounts - Normalised counts in Uncut.

Table 3 shows the total number of differentialy expressed genes between conditions under user specified cut-offs.

logFC_Down <- -1
logFC_Up <- 1
padjCutOff <- 0.05
l <- data.frame(Comparison=c(paste0("Up_In_",conditionA," vs ",conditionB),
                             paste0("Up_In_",conditionA," vs ",conditionB," (logFC>",logFC_Up,")"),
                             paste0("Up_In_",conditionB," vs ",conditionA),
                             paste0("Up_In_",conditionB," vs ",conditionA," (logFC>",logFC_Down*-1,")")),
    NumberOfGenes=c(sum(DEtable$padj < padjCutOff & !is.na(DEtable$padj) & DEtable$log2FoldChange > 0),
    sum(DEtable$padj < padjCutOff & !is.na(DEtable$padj) & DEtable$log2FoldChange > logFC_Up),
    sum(DEtable$padj < padjCutOff & !is.na(DEtable$padj) & DEtable$log2FoldChange < 0),
    sum(DEtable$padj < padjCutOff & !is.na(DEtable$padj) & DEtable$log2FoldChange < logFC_Down)
  ))

DT::datatable(l,caption = paste0("Table 3 - Number of genes passing cutoffs (padj <",padjCutOff,")"))

DEtable %>% 
   filter((padj < padjCutOff & !is.na(padj)) & ((DEtable$log2FoldChange < logFC_Down)|(DEtable$log2FoldChange > logFC_Up))) %>% 
   rio::export(file = file.path(outdir,paste0(ComparisonName,"_DGEanalysis.xlsx")))

DEtable %>% 
   filter((padj < padjCutOff & !is.na(padj)) & ((DEtable$log2FoldChange < logFC_Down)|(DEtable$log2FoldChange > logFC_Up))) %>% 
   write.table(file = file.path(outdir,paste0(ComparisonName,"_DGEanalysis.csv")),sep=",",row.names = FALSE)

DEtable %>% 
   write.table(file = file.path(outdir,paste0(ComparisonName,"_DGEanalysisUnffiltered.csv")),sep=",",row.names = FALSE)

Differential expression plots.

In this section we review two quality control plots to assess the expression level, direction and degree of change between Cut and Uncut, as well as significance of this change across all measured genes.

Volcano plot

The Volcano plot in Figure 3 shows the log2 fold change for Cut minus Uncut along the X-axis and the significance of change (the absolute Wald statistic) on the Y-axis.

ComparisonName <- paste0(conditionA,"_minus_",conditionB)    



toPlot <- c(as.vector(DEtable[order(DEtable$stat),"rowname"])[1:5],
            as.vector(DEtable[order(DEtable$stat,decreasing = TRUE),"rowname"])[1:5])

# selectedUp <- c("MEST")
# selectedDown <- c("A2M","HIST1H2BK")

DEtable$Cut_Off <- factor(ifelse((DEtable$padj < padjCutOff & !is.na(DEtable$padj)) & ((DEtable$log2FoldChange < logFC_Down)|(DEtable$log2FoldChange > logFC_Up)),"Pass",
                             "Fail"),levels = c("Fail","Pass"))


p <- ggplot(DEtable,aes(x=log2FoldChange,y=abs(stat)))+
  geom_point(aes(alpha=Cut_Off))+
  theme_bw()+
  ggtitle(ComparisonName)


p2 <- p+geom_label_repel(data = DEtable[DEtable$rowname %in% toPlot  & DEtable$padj <0.05,], 
                        aes(label=rowname), hjust=0, vjust=0, size=5)+
  xlim(c(-10,10))+
  scale_alpha_manual(values=c(Fail=0.1,Pass=1))

MA plot

The MA plot in Figure 4 shows the log2 fold change for Cut minus Uncut along the Y-axis and the log2 baseMean across all samples on the X-axis.

p <- ggplot(DEtable,aes(x=log2FoldChange,y=log2(baseMean)))+geom_point(aes(alpha=Cut_Off))+theme_bw()+ggtitle(ComparisonName)
p2 <- p+geom_label_repel(data = DEtable[DEtable$rowname %in% toPlot  & DEtable$padj <0.05,],aes(label=rowname),hjust=0,vjust=0,size=5)+xlim(c(-10,10))+scale_alpha_manual(values=c(Fail=0.1,Pass=1))+coord_flip()

Sample and Gene Heatmap

Figure 5 shows the per gene median scaled log2 expression values for significantly differential genes between Cut and Uncut.

sigGenes <- DEtable %>% 
                filter((padj < padjCutOff & !is.na(padj)) & ((DEtable$log2FoldChange < logFC_Down)|(DEtable$log2FoldChange > logFC_Up))) %>% 
                dplyr::pull(rowname)
if(length(sigGenes) > 2){
  exprsForHM <- assay(normExprs)
  exprsForHM <- exprsForHM[rownames(exprsForHM) %in% sigGenes,]
  
  
  pheatmap(t(scale(t(exprsForHM),center = TRUE,scale = FALSE)),
           clustering_distance_rows="correlation",
           clustering_distance_cols="correlation",
           annotation_col=as.data.frame(colData(normExprs))  %>% dplyr::select(!contains("sizeFactor"))  %>% dplyr::select(Group),
           scale="none",show_rownames = FALSE)
}

Figure 6 shows the per gene Z-score scaled expression values for significantly differential genes between Cut and Uncut. Z-score is the number of standard deviations from the mean, so will reveal patterns in gene expression, but not show absolute values of gene expression.

sigGenes <- DEtable %>% 
                filter((padj < padjCutOff & !is.na(padj)) & ((DEtable$log2FoldChange < logFC_Down)|(DEtable$log2FoldChange > logFC_Up))) %>% 
                dplyr::pull(rowname)
if(length(sigGenes) > 2){
  exprsForHM <- assay(normExprs)
  exprsForHM <- exprsForHM[rownames(exprsForHM) %in% sigGenes,]
  
  pheatmap(exprsForHM,
           clustering_distance_rows="correlation",
           clustering_distance_cols="correlation",
           annotation_col=as.data.frame(colData(normExprs))  %>% dplyr::select(!contains("sizeFactor"))  %>% dplyr::select(Group),
           scale="row",show_rownames = FALSE)
}

GSEA analysis of GO terms.

Following the tests of differential expression, we now perform a gene set enrichment analysis using the fgsea package to identify GO genesets which are upregulated or downregulated in comparison Cut minus Uncut.

Table 4 and figure 6 show the ranking of the differentially regulated gene sets and the gene set netowrk plots for top 40 differentially regulated by p-value.

For information on using fgsea, visit the package the package home page or our course material here

require(fgsea)
require(dplyr)
require(tibble)
require(magrittr)
require(GO.db)
require(clusterProfiler)
# GOmap <- "~/Desktop/Projects/heintz/monkeyGORefs/MfascIDToGO.txt"
# # DEtable <- read.delim("~/Desktop/Projects/brc/renvTest/AnotherTest2/DE_Genes/salmon/Group_shPTBP1_53_minus_TamDEG.xls")
# DEtable <- read.delim("/Volumes/MRCbackup/eric_mFasc/DE_Genes/Group_contra_ctx_input_minus_contra_ctx_IPDEG.xls")
# annoS <- read.delim("/Volumes/MRCbackup/eric_mFasc/MfascicularisIDtoSymbolforPipe.txt",sep="\t")

temp <- read.delim(GOsets,row.names = NULL)
# temp <- merge(temp,annoS,by=1,all=TRUE)

# temp <- temp[!temp$GOIDs == "" & !is.na(temp$Symbols),]
temp <- temp[!temp$Category == "",] %>%
  dplyr::select(TERM=Category,GENE=Genes)

# mttr <- split(as.vector(temp$Genes), as.vector(temp$Category))
# # mttr <- split(as.vector(temp$Symbols), as.vector(temp$GOIDs))
#
# mttr <- lapply(mttr,unique)
statTest <- DEtable %>% filter((!is.na(padj))) %>% dplyr::pull(stat)

# names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,2]
names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,1]

re <- as.data.frame(GO.db::GOTERM)[,-1] %>%
  dplyr::select(go_id,Term) %>%
  filter(!duplicated(go_id))

gsea_output <- GSEA(sort(statTest,decreasing = TRUE),
                    TERM2GENE = temp,
                    TERM2NAME = re,
                    seed = TRUE,
                    pvalueCutoff = 1,
                    eps = 1e-100)

gsea_output <- pairwise_termsim(gsea_output,showCategory = 40)
network_plot <- emapplot(gsea_output,cex_label_category = 0.7,showCategory = 40,color="NES") +
        ggtitle(paste0("\n",
                       "Over representation analysis for ",comparisonname))
network_plot

datatable(dplyr::select(gsea_output@result %>% arrange(pvalue),
                                     -leading_edge,
                                     -core_enrichment
  ))

gsea_output@result %>% arrange(pvalue) %>%
   rio::export(file = file.path(outdir,paste0(ComparisonName,"_GSEA.xlsx")))

GSEA analysis of KEGG terms.

Following the tests of differential expression, we now perform a gene set enrichment analysis using the fgsea package to identify GO genesets which are upregulated or downregulated in comparison Cut minus Uncut.

Table 5 and figure 7 show the ranking of the differentially regulated gene sets and the gene set netowrk plots for top 40 differentially regulated by p-value.

For information on using fgsea, visit the package the package home page or our course material here

require(fgsea)
require(dplyr)
require(tibble)
require(magrittr)
require(GO.db)
require(clusterProfiler)
# GOmap <- "~/Desktop/Projects/heintz/monkeyGORefs/MfascIDToGO.txt"
# # DEtable <- read.delim("~/Desktop/Projects/brc/renvTest/AnotherTest2/DE_Genes/salmon/Group_shPTBP1_53_minus_TamDEG.xls")
# DEtable <- read.delim("/Volumes/MRCbackup/eric_mFasc/DE_Genes/Group_contra_ctx_input_minus_contra_ctx_IPDEG.xls")
# annoS <- read.delim("/Volumes/MRCbackup/eric_mFasc/MfascicularisIDtoSymbolforPipe.txt",sep="\t")

temp <- read.delim(KEGGsets,row.names = NULL)
# temp <- merge(temp,annoS,by=1,all=TRUE)

# temp <- temp[!temp$GOIDs == "" & !is.na(temp$Symbols),]
temp <- temp[!temp$Category == "",] %>% 
  dplyr::select(TERM=Category,GENE=Genes)

# mttr <- split(as.vector(temp$Genes), as.vector(temp$Category))
# # mttr <- split(as.vector(temp$Symbols), as.vector(temp$GOIDs))
# 
# mttr <- lapply(mttr,unique)
statTest <- DEtable %>% filter((!is.na(padj))) %>% dplyr::pull(stat)

# names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,2]
names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,1]

re <- as.data.frame(KEGGREST::keggList("pathway")) %>% rownames_to_column("kegg_id") %>% 
  set_colnames(c("KEGGID","TERM")) %>% 
  filter(!duplicated(KEGGID))

temp <- temp[temp[,1] %in% re[,1],]
gsea_output <- GSEA(sort(statTest,decreasing = TRUE),
                    TERM2GENE = temp,
                    TERM2NAME = re,
                    seed = TRUE,
                    pvalueCutoff = 1,
                    eps = 1e-100)
    
gsea_output <- pairwise_termsim(gsea_output,showCategory = 40)
network_plot <- emapplot(gsea_output,cex_label_category = 0.7,showCategory = 40,color="NES") +
        ggtitle(paste0("\n",
                       "Over representation analysis for ",comparisonname))
network_plot

datatable(dplyr::select(gsea_output@result %>% arrange(pvalue),
                                     -leading_edge,
                                     -core_enrichment
  ))

gsea_output@result %>% arrange(pvalue) %>% 
   rio::export(file = file.path(outdir,paste0(ComparisonName,"_KEGG_GSEA.xlsx")))

GSEA analysis of PFAM terms.

Following the tests of differential expression, we now perform a gene set enrichment analysis using the fgsea package to identify GO genesets which are upregulated or downregulated in comparison Cut minus Uncut.

Table 6 and figure 8 show the ranking of the differentially regulated gene sets and the gene set netowrk plots for top 40 differentially regulated by p-value.

For information on using fgsea, visit the package the package home page or our course material here

require(fgsea)
require(dplyr)
require(tibble)
require(magrittr)
require(GO.db)
require(clusterProfiler)
# GOmap <- "~/Desktop/Projects/heintz/monkeyGORefs/MfascIDToGO.txt"
# # DEtable <- read.delim("~/Desktop/Projects/brc/renvTest/AnotherTest2/DE_Genes/salmon/Group_shPTBP1_53_minus_TamDEG.xls")
# DEtable <- read.delim("/Volumes/MRCbackup/eric_mFasc/DE_Genes/Group_contra_ctx_input_minus_contra_ctx_IPDEG.xls")
# annoS <- read.delim("/Volumes/MRCbackup/eric_mFasc/MfascicularisIDtoSymbolforPipe.txt",sep="\t")

temp <- read.delim(PFAMsets,row.names = NULL)
# temp <- merge(temp,annoS,by=1,all=TRUE)

# temp <- temp[!temp$GOIDs == "" & !is.na(temp$Symbols),]
temp <- temp[!temp$Category == "",] %>% 
  dplyr::select(TERM=Category,GENE=Genes)

# mttr <- split(as.vector(temp$Genes), as.vector(temp$Category))
# # mttr <- split(as.vector(temp$Symbols), as.vector(temp$GOIDs))
# 
# mttr <- lapply(mttr,unique)
statTest <- DEtable %>% filter((!is.na(padj))) %>% dplyr::pull(stat)

# names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,2]
names(statTest) <- DEtable %>% filter((!is.na(padj))) %>% .[,1]

# re <- as.data.frame(KEGGREST::keggList("pathway")) %>% rownames_to_column("kegg_id") %>% 
#   set_colnames(c("KEGGID","TERM")) %>% 
#   filter(!duplicated(KEGGID))

gsea_output <- GSEA(sort(statTest,decreasing = TRUE),
                    TERM2GENE = temp,
                    seed = TRUE,
                    pvalueCutoff = 1,
                    eps = 1e-100)
    
gsea_output <- pairwise_termsim(gsea_output,showCategory = 40)
network_plot <- emapplot(gsea_output,cex_label_category = 0.7,showCategory = 40,color="NES") +
        ggtitle(paste0("\n",
                       "Over representation analysis for ",comparisonname))
network_plot

datatable(dplyr::select(gsea_output@result %>% arrange(pvalue),
                                     -leading_edge,
                                     -core_enrichment
  ))

gsea_output@result %>% arrange(pvalue) %>% 
   rio::export(file = file.path(outdir,paste0(ComparisonName,"_PFAM_GSEA.xlsx")))

Glossary of Terms

Return to table of contents

Key terms with Course Links

Software used

Return to table of contents

Packages outside of R

Package	Citation
Picard Tools	https://broadinstitute.github.io/picard/

R packages

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## Random number generation:
##  RNG:     L'Ecuyer-CMRG 
##  Normal:  Inversion 
##  Sample:  Rejection 
##  
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] Rsamtools_2.10.0             Biostrings_2.62.0           
##  [3] XVector_0.34.0               IsoformSwitchAnalyzeR_1.16.0
##  [5] DEXSeq_1.40.0                limma_3.50.1                
##  [7] GenomicFeatures_1.46.5       org.Dr.eg.db_3.14.0         
##  [9] org.Dm.eg.db_3.14.0          BiocParallel_1.28.3         
## [11] aws.s3_0.3.21                futile.logger_1.4.3         
## [13] gridExtra_2.3                GO.db_3.14.0                
## [15] magick_2.7.3                 org.Mm.eg.db_3.14.0         
## [17] org.Hs.eg.db_3.14.0          enrichplot_1.14.2           
## [19] clusterProfiler_4.2.2        msigdbr_7.5.1               
## [21] ComplexHeatmap_2.10.0        highcharter_0.9.4           
## [23] GlossyBRC_0.1.0              jsonlite_1.8.0              
## [25] rjson_0.2.21                 rio_0.5.29                  
## [27] plotly_4.10.0                RColorBrewer_1.1-3          
## [29] htmltools_0.5.2              knitr_1.38                  
## [31] rmarkdown_2.13               kableExtra_1.3.4            
## [33] ggrepel_0.9.1                DT_0.22                     
## [35] fgsea_1.20.0                 goseq_1.46.0                
## [37] geneLenDataBase_1.30.0       BiasedUrn_1.07              
## [39] GSEABase_1.56.0              graph_1.72.0                
## [41] annotate_1.72.0              XML_3.99-0.8                
## [43] AnnotationDbi_1.56.2         pheatmap_1.0.12             
## [45] tximport_1.22.0              crosstalk_1.2.0             
## [47] DESeq2_1.34.0                SummarizedExperiment_1.24.0 
## [49] Biobase_2.54.0               MatrixGenerics_1.6.0        
## [51] matrixStats_0.61.0           GenomicRanges_1.46.1        
## [53] GenomeInfoDb_1.30.1          IRanges_2.28.0              
## [55] S4Vectors_0.32.3             BiocGenerics_0.40.0         
## [57] forcats_0.5.1                stringr_1.4.0               
## [59] dplyr_1.0.8                  purrr_0.3.4                 
## [61] readr_2.1.2                  tidyr_1.2.0                 
## [63] tibble_3.1.7                 ggplot2_3.3.5               
## [65] tidyverse_1.3.1              magrittr_2.0.2              
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3                rtracklayer_1.54.0           
##   [3] bit64_4.0.5                   DelayedArray_0.20.0          
##   [5] hwriter_1.3.2                 data.table_1.14.2            
##   [7] AnnotationFilter_1.18.0       KEGGREST_1.34.0              
##   [9] RCurl_1.98-1.6                doParallel_1.0.17            
##  [11] generics_0.1.2                lambda.r_1.2.4               
##  [13] RSQLite_2.2.10                shadowtext_0.1.2             
##  [15] bit_4.0.4                     tzdb_0.2.0                   
##  [17] rlist_0.4.6.2                 httpuv_1.6.5                 
##  [19] webshot_0.5.2                 xml2_1.3.3                   
##  [21] lubridate_1.8.0               assertthat_0.2.1             
##  [23] viridis_0.6.2                 xfun_0.39                    
##  [25] jquerylib_0.1.4               hms_1.1.1                    
##  [27] promises_1.2.0.1              babelgene_22.3               
##  [29] evaluate_0.15                 fansi_1.0.2                  
##  [31] restfulr_0.0.13               progress_1.2.2               
##  [33] dbplyr_2.1.1                  readxl_1.3.1                 
##  [35] igraph_1.3.0                  DBI_1.1.2                    
##  [37] geneplotter_1.72.0            quantmod_0.4.20              
##  [39] htmlwidgets_1.5.4             ellipsis_0.3.2               
##  [41] ggnewscale_0.4.7              backports_1.4.1              
##  [43] biomaRt_2.50.3                vctrs_0.4.1                  
##  [45] ensembldb_2.18.4              TTR_0.24.3                   
##  [47] cachem_1.0.6                  withr_2.5.0                  
##  [49] aws.signature_0.6.0           ggforce_0.3.3                
##  [51] BSgenome_1.62.0               vroom_1.5.7                  
##  [53] treeio_1.18.1                 GenomicAlignments_1.30.0     
##  [55] xts_0.12.2                    prettyunits_1.1.1            
##  [57] svglite_2.1.0                 cluster_2.1.2                
##  [59] DOSE_3.20.1                   ape_5.6-2                    
##  [61] lazyeval_0.2.2                crayon_1.5.0                 
##  [63] genefilter_1.76.0             labeling_0.4.2               
##  [65] edgeR_3.36.0                  pkgconfig_2.0.3              
##  [67] tweenr_1.0.2                  ProtGenerics_1.26.0          
##  [69] nlme_3.1-153                  pals_1.7                     
##  [71] rlang_1.0.6                   lifecycle_1.0.1              
##  [73] downloader_0.4                filelock_1.0.2               
##  [75] BiocFileCache_2.2.1           modelr_0.1.8                 
##  [77] AnnotationHub_3.2.2           dichromat_2.0-0.1            
##  [79] VennDiagram_1.7.3             cellranger_1.1.0             
##  [81] polyclip_1.10-0               Matrix_1.3-4                 
##  [83] aplot_0.1.6                   tximeta_1.12.4               
##  [85] zoo_1.8-9                     base64enc_0.1-3              
##  [87] reprex_2.0.1                  GlobalOptions_0.1.2          
##  [89] png_0.1-7                     viridisLite_0.4.0            
##  [91] bitops_1.0-7                  blob_1.2.2                   
##  [93] shape_1.4.6                   qvalue_2.26.0                
##  [95] gridGraphics_0.5-1            scales_1.1.1                 
##  [97] memoise_2.0.1                 plyr_1.8.7                   
##  [99] zlibbioc_1.40.0               scatterpie_0.1.7             
## [101] compiler_4.1.2                BiocIO_1.4.0                 
## [103] clue_0.3-61                   cli_3.6.0                    
## [105] patchwork_1.1.1               formatR_1.12                 
## [107] MASS_7.3-54                   mgcv_1.8-38                  
## [109] tidyselect_1.1.2              stringi_1.7.6                
## [111] highr_0.9                     yaml_2.3.5                   
## [113] GOSemSim_2.20.0               locfit_1.5-9.4               
## [115] sass_0.4.1                    fastmatch_1.1-3              
## [117] tools_4.1.2                   parallel_4.1.2               
## [119] circlize_0.4.15               rstudioapi_0.13              
## [121] foreach_1.5.2                 foreign_0.8-81               
## [123] farver_2.1.0                  ggraph_2.0.5                 
## [125] BiocManager_1.30.16           digest_0.6.29                
## [127] shiny_1.7.2                   Rcpp_1.0.8                   
## [129] broom_0.7.12                  later_1.3.0                  
## [131] BiocVersion_3.14.0            httr_1.4.2                   
## [133] colorspace_2.0-3              rvest_1.0.2                  
## [135] fs_1.5.2                      splines_4.1.2                
## [137] statmod_1.4.36                yulab.utils_0.0.4            
## [139] tidytree_0.3.9                graphlayouts_0.8.0           
## [141] mapproj_1.2.9                 ggplotify_0.1.0              
## [143] systemfonts_1.0.4             xtable_1.8-4                 
## [145] futile.options_1.0.1          ggtree_3.2.1                 
## [147] tidygraph_1.2.0               ggfun_0.0.6                  
## [149] R6_2.5.1                      DRIMSeq_1.22.0               
## [151] mime_0.12                     pillar_1.7.0                 
## [153] glue_1.6.2                    fastmap_1.1.0                
## [155] interactiveDisplayBase_1.32.0 codetools_0.2-18             
## [157] maps_3.4.1                    utf8_1.2.2                   
## [159] bslib_0.3.1                   lattice_0.20-45              
## [161] curl_4.3.2                    zip_2.2.0                    
## [163] openxlsx_4.2.5                survival_3.2-13              
## [165] munsell_0.5.0                 DO.db_2.9                    
## [167] GetoptLong_1.0.5              GenomeInfoDbData_1.2.7       
## [169] iterators_1.0.14              haven_2.4.3                  
## [171] reshape2_1.4.4                gtable_0.3.0